const Array2DRowRealMatrix = require('./Array2DRowRealMatrix'),
  LUDecomposition = require('./LUDecomposition');

function FastHessian(integralImage, octaves, initSample, threshold, balanceValue) {
  this.mIntegralImage = integralImage;
  this.mOctaves = octaves;
  this.mInitSample = initSample;
  this.mThreshold = threshold;
  this.mWidth = integralImage.getWidth();
  this.mHeight = integralImage.getHeight();
  this.mInterestPoints = null;
  this.mRecalculateInterestPoints = true;
  this.mLayers = null;

  var filter_map = [[
    0,
    1,
    2,
    3
  ],
    [
      1,
      3,
      4,
      5
    ],
    [
      3,
      5,
      6,
      7
    ],
    [
      5,
      7,
      8,
      9
    ],
    [
      7,
      9,
      10,
      11
    ]
  ];

  this.getIPoints = function () {
    if (this.mInterestPoints == null || this.mRecalculateInterestPoints) {
      this.mInterestPoints = [];
      this.buildResponseMap();
      var b,
        m,
        t;
      for (var o = 0; o < this.mOctaves; o++) {
        for (var i = 0; i <= 1; i++) {
          b = this.mLayers[filter_map[o][i]];
          m = this.mLayers[filter_map[o][i + 1]];
          t = this.mLayers[filter_map[o][i + 2]];

          // loop over middle response layer at density of the most
          // sparse layer (always top), to find maxima across scale and space
          for (var r = 0; r < t.getHeight(); r++) {
            for (var c = 0; c < t.getWidth(); c++) {
              if (this.isExtremum(r, c, t, m, b)) {
                var point = this.interpolateExtremum(r, c, t, m, b);
                if (point != null) {
                  this.mInterestPoints.push(point);
                }
              }
            }
          }
        }
      }
    }
    return this.mInterestPoints;
  }

  this.buildResponseMap = function () {
    this.mLayers = [];
    var w = Math.floor(this.mWidth / this.mInitSample);
    var h = Math.floor(this.mHeight / this.mInitSample);
    var s = this.mInitSample;
    const ResponseLayer = require('./ResponseLayer');
    if (this.mOctaves >= 1) {
      this.mLayers.push(new ResponseLayer(w, h, s, 9, this.mIntegralImage));
      this.mLayers.push(new ResponseLayer(w, h, s, 15, this.mIntegralImage));
      this.mLayers.push(new ResponseLayer(w, h, s, 21, this.mIntegralImage));
      this.mLayers.push(new ResponseLayer(w, h, s, 27, this.mIntegralImage));
    }
    if (this.mOctaves >= 2) {
      this.mLayers.push(new ResponseLayer(Math.floor(w / 2), Math.floor(h / 2), s * 2, 39, this.mIntegralImage));
      this.mLayers.push(new ResponseLayer(Math.floor(w / 2), Math.floor(h / 2), s * 2, 51, this.mIntegralImage));
    }
    if (this.mOctaves >= 3) {
      this.mLayers.push(new ResponseLayer(Math.floor(w / 4), Math.floor(h / 4), s * 4, 75, this.mIntegralImage));
      this.mLayers.push(new ResponseLayer(Math.floor(w / 4), Math.floor(h / 4), s * 4, 99, this.mIntegralImage));
    }
    if (this.mOctaves >= 4) {
      this.mLayers.push(new ResponseLayer(Math.floor(w / 8), Math.floor(h / 8), s * 8, 147, this.mIntegralImage));
      this.mLayers.push(new ResponseLayer(Math.floor(w / 8), Math.floor(h / 8), s * 8, 195, this.mIntegralImage));
    }
    if (this.mOctaves >= 5) {
      this.mLayers.push(new ResponseLayer(Math.floor(w / 16), Math.floor(h / 16), s * 16, 291, this.mIntegralImage));
      this.mLayers.push(new ResponseLayer(Math.floor(w / 16), Math.floor(h / 16), s * 16, 387, this.mIntegralImage));
    }
  }

  this.isExtremum = function (r, c, t, m, b) {
    var layerBorder = Math.floor((t.getFilter() + 1) / (2 * t.getStep()));
    if (r <= layerBorder || r >= t.getHeight() - layerBorder || c <= layerBorder || c >= t.getWidth() - layerBorder) {
      return false;
    }
    var candidate = m.getResponse(r, c, t);
    if (candidate < this.mThreshold)
      return false;

    //See if the response in 3x3x3 is greater, then it isn't a local maxima
    for (var rr = -1; rr <= 1; rr++) {
      for (var cc = -1; cc <= 1; cc++) {
        if (r + rr < 0 || r + rr >= t.getHeight() || c + cc < 0 || c + cc >= t.getWidth())
          return false;
        if (t.getResponse(r + rr, c + cc) >= candidate || ((rr != 0 || cc != 0) && m.getResponse(r + rr, c + cc, t) >= candidate) || b.getResponse(r + rr, c + cc, t) >= candidate)
          return false;
      }
    }
    return true;
  }

  this.interpolateExtremum = function (r, c, t, m, b) {
    //should check to make sure that m's filter value is less than t's and greater than b's
    var filterStep = m.getFilter() - b.getFilter();
    var xi = 0,
      xr = 0,
      xc = 0;
    var values = this.interpolateStep(r, c, t, m, b);
    xi = values[0];
    xr = values[1];
    xc = values[2];
    if (Math.abs(xi) < 0.5 && Math.abs(xr) < 0.5 && Math.abs(xc) < 0.5) {
      //WE"VE GOT AN INTEREST POINT HERE
      var x = (c + xc) * t.getStep();
      var y = (r + xr) * t.getStep();
      var scale = 0.1333 * (m.getFilter() + xi * filterStep);
      var laplacian = Math.floor(m.getLaplacian(r, c, t));
      const SURFInterestPoint = require('./SURFInterestPoint');
      return new SURFInterestPoint(x, y, scale, laplacian);
    }
    return null;
  }

  this.interpolateStep = function (r, c, t, m, b) {
    var values = new Array(3);
    var partialDerivs = this.getPartialDerivativeMatrix(r, c, t, m, b);
    var hessian3D = this.getHessian3DMatrix(r, c, t, m, b);
    var solver = new LUDecomposition(hessian3D).getSolver();
    var X = solver.getInverse().multiply(partialDerivs);

    //values of them are used
    //xi
    values[0] = -X.getEntry(2, 0);
    //xr
    values[1] = -X.getEntry(1, 0);
    //xc
    values[2] = -X.getEntry(0, 0);

    return values;
  }

  this.getPartialDerivativeMatrix = function (r, c, t, m, b) {
    var derivs = new Array(3);
    derivs[0] = new Array(1);
    derivs[1] = new Array(1);
    derivs[2] = new Array(1);
    derivs[0][0] = (m.getResponse(r, c + 1, t) - m.getResponse(r, c - 1, t)) / 2.0;
    derivs[1][0] = (m.getResponse(r + 1, c, t) - m.getResponse(r - 1, c, t)) / 2.0;
    derivs[2][0] = (t.getResponse(r, c) - b.getResponse(r, c, t)) / 2.0;
    var matrix = new Array2DRowRealMatrix(derivs);
    return matrix;
  }

  this.getHessian3DMatrix = function (r, c, t, m, b) {
    //Layout:
    //  [dxx][dxy][dxs]
    //  [dxy][dyy][dys]
    //  [dxs][dys][dss]
    var hessian = new Array(3);
    for (var i = 0; i < hessian.length; i++) {
      hessian[i] = new Array(3)
    }

    var v = m.getResponse(r, c, t);

    //dxx
    hessian[0][0] = m.getResponse(r, c + 1, t) +
      m.getResponse(r, c - 1, t) -
      2 * v;

    //dyy
    hessian[1][1] = m.getResponse(r + 1, c, t) +
      m.getResponse(r - 1, c, t) -
      2 * v;

    //dss
    hessian[2][2] = t.getResponse(r, c) +
      b.getResponse(r, c, t) -
      2 * v;

    //dxy
    hessian[0][1] = hessian[1][0] = (m.getResponse(r + 1, c + 1, t) -
      m.getResponse(r + 1, c - 1, t) -
      m.getResponse(r - 1, c + 1, t) +
      m.getResponse(r - 1, c - 1, t)) / 4.0;

    //dxs
    hessian[0][2] = hessian[2][0] = (t.getResponse(r, c + 1) -
      t.getResponse(r, c - 1) -
      b.getResponse(r, c + 1, t) +
      b.getResponse(r, c - 1, t)) / 4.0;

    //dys
    hessian[1][2] = hessian[2][1] = (t.getResponse(r + 1, c) -
      t.getResponse(r - 1, c) -
      b.getResponse(r + 1, c, t) +
      b.getResponse(r - 1, c, t)) / 4.0;

    return new Array2DRowRealMatrix(hessian);
  }
}

module.exports = FastHessian
